##################################################
#
# COMPARISONS OF MODELS BY RELIGIOUS PRESENCE IN
# WINNING COALITION
#
# Paul Musgrave and Yu-Ming Liou
# musgrave@umass.edu and yl254@georgetown.edu
#
# Created 27 February 2016
#
##################################################

##################################################
# R Housekeeping
##################################################

# Remove and empty workspace
rm(list=ls())

#  #  #  #  #  #  #  #  #  #  #  #  #  #  #  #  # 
# Load Packages
#  #  #  #  #  #  #  #  #  #  #  #  #  #  #  #  # 

library(foreign)
library(sandwich)
library(apsrtable)
library(tonymisc)  # includes summaryR()
library(MASS) # for polr() which runs our ologit

#  #  #  #  #  #  #  #  #  #  #  #  #  #  #  #  # 
# Load Functions
#  #  #  #  #  #  #  #  #  #  #  #  #  #  #  #  # 

# Change this to  your working directory
mywd <- "/Users/paulmusgrave/Dropbox/0001 Academic Projects/Ongoing/0127 Oil Islam Women/ISQ Accepted Submission/Replication"

setwd(paste(mywd,"/Code",sep=""))
source("statifyFN.R")
source("multipleOLSFN.R")
source("multipleOlogitFN.R")
source("genModelsFN.R")
source("plotlineFN.R")
source("plotprepFN.R")
source("modelDefinitions.R")


#  #  #  #  #  #  #  #  #  #  #  #  #  #  #  #  # 
# Load Data
#  #  #  #  #  #  #  #  #  #  #  #  #  #  #  #  # 

setwd(paste(mywd,"/Data",sep=""))

data <- read.dta("NewDVModelsXSforISQFINAL2015.dta")

data$oil_gas_pc_11_log <- log(data$oil_gas_pc_11+1)
data$oil_gas_pc_07_log <- log(data$oil_gas_pc_07+1)
data$oil_gas_pc_04_log <- log(data$oil_gas_pc_04+1)
data$oil_gas_pc_02_log <- log(data$oil_gas_pc_02+1)
data$oil_gas_pc_09_log <- log(data$oil_gas_pc_09+1)


data$oil_gas_pc_11 <- data$oil_gas_pc_11 / 1000
data$oil_gas_pc_07 <- data$oil_gas_pc_07 / 1000
data$oil_gas_pc_04 <- data$oil_gas_pc_04 / 1000
data$oil_gas_pc_02 <- data$oil_gas_pc_02 / 1000
data$oil_gas_pc_09 <- data$oil_gas_pc_09 / 1000


#################################################
# Regressions
#################################################

### ALCOHOL ###
out.alcohol.1 <- multipleOLS(formulas=genModels(c("totalalcohollitres"),
                                             c(p140.04,p142.04)),
                                             y=data[data$taxsupport_05 ==1 &
                                             		data$p_polity2_04 <= 5,],
                          modelnames=c(names(p140.04),names(p142.04)))

out.alcohol.0 <- multipleOLS(formulas=genModels(c("totalalcohollitres"),
                                             c(p140.04,p142.04)),
                                             y=data[data$taxsupport_05 ==0 &
                                             		data$p_polity2_04 <= 5,],
                          modelnames=c(names(p140.04),names(p142.04)))
                          
                          
### MALE SMOKING ###

out.malesmoking.1 <- multipleOLS(formulas=genModels(c("malesmoking"),
                                            c(p142.09)),
                                             y=data[data$taxsupport_08 ==1 &
                                             		data$p_polity2_09 <= 5,],
				          	modelnames=c(names(p142.09)))

out.malesmoking.0 <- multipleOLS(formulas=genModels(c("malesmoking"),
                                            c(p142.09)),
                                             y=data[data$taxsupport_08 ==0 &
                                             		data$p_polity2_09 <= 5,],
			
				          	modelnames=c(names(p142.09)))

### FEMALE SMOKING ###
out.femsmoking.1 <- multipleOLS(formulas=genModels(c("femsmoking"),
                                             c(p142.09)),
                                             y=data[data$taxsupport_08 ==1 &
                                             		data$p_polity2_09 <= 5,],
			
                          modelnames=c(names(p142.09)))
out.femsmoking.0 <- multipleOLS(formulas=genModels(c("femsmoking"),
                                             c(p142.09)),
                                             y=data[data$taxsupport_08 ==0 &
                                             		data$p_polity2_09 <= 5,],
			
                          modelnames=c(names(p142.09)))
                          

### Total Fertility ###
out.totfert.1 <- multipleOLS(formulas=genModels(c("totfert"),
                                            c(p140alt,p142base)),
                                             y=data[data$taxsupport_01 ==1 &
                                      data$p_polity2_02 <= 5,],
				          	modelnames=c(names(p140alt),names(p142base)))  
				          	                                   
out.totfert.0 <- multipleOLS(formulas=genModels(c("totfert"),
                                            c(p140alt,p142base)),
                                             y=data[data$taxsupport_01 ==0 &
                                      data$p_polity2_02 <= 5,],
				          	modelnames=c(names(p140alt),names(p142base)))  
				          	

### Adolescent Fertility ###
                                    
out.adfert.1 <- multipleOLS(formulas=genModels(c("adfert"),
                                            c(p140alt,p142base)),
                                             y=data[data$taxsupport_01 ==1 &
                                      data$p_polity2_02 <= 5,],
				          	modelnames=c(names(p140alt),names(p142base)))  
				          	                                   
out.adfert.0 <- multipleOLS(formulas=genModels(c("adfert"),
                                            c(p140alt,p142base)),
                                             y=data[data$taxsupport_01 ==0 &
                                      data$p_polity2_02 <= 5,],
				          	modelnames=c(names(p140alt),names(p142base)))                                                     
                                                                                                                                
 
### GII 11 ###

out.GII.0 <- multipleOLS(formulas=genModels(c("GII_UN_11"),
                                                 c(p140.11,
                                                   p142.11)),
                              modelnames=c(names(p140.11),
                                           names(p142.11)),
                              y=data[data$taxsupport_08 ==0 &
                                      data$p_polity2_11 <= 5,])

out.GII.1 <- multipleOLS(formulas=genModels(c("GII_UN_11"),
                                                 c(p140.11,
                                                   p142.11)),
                              modelnames=c(names(p140.11),
                                           names(p142.11)),
                              y=data[data$taxsupport_08 ==1 &
                                      data$p_polity2_11 <= 5,])    



### Contraceptive Prevalence ###

out.contraception.1 <- multipleOLS(formulas=genModels(c("avgcontraception"),
                                               c(p140alt,p142base)),
                            y=data[data$taxsupport_01 ==1 &
                                     data$p_polity2_02 <= 5,],
                            modelnames=c(names(p140alt),names(p142base)))  

out.contraception.0 <- multipleOLS(formulas=genModels(c("avgcontraception"),
                                               c(p140alt,p142base)),
                            y=data[data$taxsupport_01 ==0 &
                                     data$p_polity2_02 <= 5,],
                            modelnames=c(names(p140alt),names(p142base)))                                                     


#################################################
# Plotting vars
#################################################

vars.to.plot <- c("oil_gas_pc_11","oil_gas_pc_09","oil_gas_pc_07",
                  "oil_gas_pc_04","oil_gas_pc_02","oil_gas_valuePOPred")


# Adfert
plot.adfert.1 <- do.call(rbind.data.frame,lapply(out.adfert.1,
                                                 plotprep,
                                                 vars.to.plot))
plot.adfert.1 <- data.frame(plot.adfert.1,c("Religious"))
colnames(plot.adfert.1)[ncol(plot.adfert.1)] <- c("Type")
plot.adfert.1$Offset <- -.1


plot.adfert.0 <- do.call(rbind.data.frame,lapply(out.adfert.0,
                                                 plotprep,
                                                 vars.to.plot))
plot.adfert.0 <- data.frame(plot.adfert.0,c("Secular"))
colnames(plot.adfert.0)[ncol(plot.adfert.0)] <- c("Type")
plot.adfert.0$Offset <- .1


# totfert
plot.totfert.1 <- do.call(rbind.data.frame,lapply(out.totfert.1,
                                                  plotprep,
                                                  vars.to.plot))
plot.totfert.1 <- data.frame(plot.totfert.1,c("Religious"))
colnames(plot.totfert.1)[ncol(plot.totfert.1)] <- c("Type")
plot.totfert.1$Offset <- -.1


plot.totfert.0 <- do.call(rbind.data.frame,lapply(out.totfert.0,
                                                  plotprep,
                                                  vars.to.plot))
plot.totfert.0 <- data.frame(plot.totfert.0,c("Secular"))
colnames(plot.totfert.0)[ncol(plot.totfert.0)] <- c("Type")
plot.totfert.0$Offset <- .1

# GII
plot.GII.1 <- do.call(rbind.data.frame,lapply(out.GII.1,
                                              plotprep,
                                              vars.to.plot))
plot.GII.1 <- data.frame(plot.GII.1,c("Religious"))
colnames(plot.GII.1)[ncol(plot.GII.1)] <- c("Type")
plot.GII.1$Offset <- -.1


plot.GII.0 <- do.call(rbind.data.frame,lapply(out.GII.0,
                                              plotprep,
                                              vars.to.plot))
plot.GII.0 <- data.frame(plot.GII.0,c("Secular"))
colnames(plot.GII.0)[ncol(plot.GII.0)] <- c("Type")
plot.GII.0$Offset <- .1




# malesmoking
plot.malesmoking.1 <- do.call(rbind.data.frame,lapply(out.malesmoking.1,
                                                      plotprep,
                                                      vars.to.plot))
plot.malesmoking.1 <- data.frame(plot.malesmoking.1,c("Religious"))
colnames(plot.malesmoking.1)[ncol(plot.malesmoking.1)] <- c("Type")
plot.malesmoking.1$Offset <- -.1


plot.malesmoking.0 <- do.call(rbind.data.frame,lapply(out.malesmoking.0,
                                                      plotprep,
                                                      vars.to.plot))
plot.malesmoking.0 <- data.frame(plot.malesmoking.0,c("Secular"))
colnames(plot.malesmoking.0)[ncol(plot.malesmoking.0)] <- c("Type")
plot.malesmoking.0$Offset <- .1



# femsmoking
plot.femsmoking.1 <- do.call(rbind.data.frame,lapply(out.femsmoking.1,
                                                     plotprep,
                                                     vars.to.plot))
plot.femsmoking.1 <- data.frame(plot.femsmoking.1,c("Religious"))
colnames(plot.femsmoking.1)[ncol(plot.femsmoking.1)] <- c("Type")
plot.femsmoking.1$Offset <- -.1


plot.femsmoking.0 <- do.call(rbind.data.frame,lapply(out.femsmoking.0,
                                                     plotprep,
                                                     vars.to.plot))
plot.femsmoking.0 <- data.frame(plot.femsmoking.0,c("Secular"))
colnames(plot.femsmoking.0)[ncol(plot.femsmoking.0)] <- c("Type")
plot.femsmoking.0$Offset <- .1



# alcohol
plot.alcohol.1 <- do.call(rbind.data.frame,lapply(out.alcohol.1,
                                                  plotprep,
                                                  vars.to.plot))
plot.alcohol.1 <- data.frame(plot.alcohol.1,c("Religious"))
colnames(plot.alcohol.1)[ncol(plot.alcohol.1)] <- c("Type")
plot.alcohol.1$Offset <- -.1


plot.alcohol.0 <- do.call(rbind.data.frame,lapply(out.alcohol.0,
                                                  plotprep,
                                                  vars.to.plot))
plot.alcohol.0 <- data.frame(plot.alcohol.0,c("Secular"))
colnames(plot.alcohol.0)[ncol(plot.alcohol.0)] <- c("Type")
plot.alcohol.0$Offset <- .1



# contraception
plot.contraception.1 <- do.call(rbind.data.frame,lapply(out.contraception.1,
                                                        plotprep,
                                                        vars.to.plot))
plot.contraception.1 <- data.frame(plot.contraception.1,c("Religious"))
colnames(plot.contraception.1)[ncol(plot.contraception.1)] <- c("Type")
plot.contraception.1$Offset <- -.1


plot.contraception.0 <- do.call(rbind.data.frame,lapply(out.contraception.0,
                                                        plotprep,
                                                        vars.to.plot))
plot.contraception.0 <- data.frame(plot.contraception.0,c("Secular"))
colnames(plot.contraception.0)[ncol(plot.contraception.0)] <- c("Type")
plot.contraception.0$Offset <- .1


## Create plotting matrices

plot.contraception <- rbind(plot.contraception.1,plot.contraception.0)
plot.alcohol <- rbind(plot.alcohol.1,plot.alcohol.0)
plot.femsmoking <- rbind(plot.femsmoking.1,plot.femsmoking.0)
plot.malesmoking <- rbind(plot.malesmoking.1,plot.malesmoking.0)
plot.GII <- rbind(plot.GII.1,plot.GII.0)
plot.totfert <- rbind(plot.totfert.1,plot.totfert.0)
plot.adfert <- rbind(plot.adfert.1,plot.adfert.0)

tuneup <- function(x){
  x$Index <- 1:(nrow(x)/2)
  x$pch <- 21
  x[x$Type == "Secular",]$pch<- 23
  x$X <- x$Index + x$Offset
  x$col <- "black"
  x$significance <- abs(x$Estimate / x$SE)
  x$bg <- "gray"
  x$bg <- ifelse (x$significance<=1.65, "white",
                  ifelse(x$significance>=1.96,"black",
                         "gray"))
  return(x)
}

plot.contraception <- tuneup(plot.contraception)
plot.alcohol <- tuneup(plot.alcohol)
plot.femsmoking <- tuneup(plot.femsmoking)
plot.malesmoking <- tuneup(plot.malesmoking)
plot.GII <- tuneup(plot.GII)
plot.totfert <- tuneup(plot.totfert)
plot.adfert <- tuneup(plot.adfert)



#  #  #  #  #  #  #  #  #  #  #  #  #  #  #  #  # 
# Graphing OIPC
#  #  #  #  #  #  #  #  #  #  #  #  #  #  #  #  # 

setwd("../Drafts/Charts")
pdf(file="ISQFINALAppendixFigure01.pdf",height=8,width=7)

par(mfrow=c(7,1),oma=c(0,0,3,1),mar=c(2,2,1,1))

plot(1,
     pch="",
     ylim=c(1.15*min(plot.totfert$Y0),1.15*max(plot.totfert$Y1)), 
     xlim=c(.5,max(plot.totfert$Index)+.5), 
     xaxt="n",
     xlab="",
     ylab="Coefficient Estimate",
     main="Total Fertility",
     cex.main=1.25,
     adj=1)
abline(h=0,lty="dotdash")
abline(v=3.5,lty="solid")

segments(y0=plot.totfert$Y0,
         y1=plot.totfert$Y1,
         x0=plot.totfert$X,
         x1=plot.totfert$X,
         col=plot.totfert$col,
         lwd=1,
         lty=1)
points(x=plot.totfert$X,
       y=plot.totfert$Estimate,
       col=plot.totfert$col,
       pch=plot.totfert$pch,
       bg=plot.totfert$bg,
       cex=1.5)
axis(side=1,at=plot.totfert[1:8,]$Index,
     labels=plot.totfert[1:8,]$Model,
     cex.axis=1.25, tick=TRUE)



plot(1,
     pch="",
     ylim=c(1.15*min(plot.adfert$Y0),1.15*max(plot.adfert$Y1)), 
     xlim=c(.5,max(plot.adfert$Index)+.5), 
     xaxt="n",
     xlab="",
     ylab="Coefficient Estimate",
     main="Adolescent Fertility",
     cex.main=1.25,
     adj=1)
abline(h=0,lty="dotdash")
abline(v=3.5,lty="solid")

segments(y0=plot.adfert$Y0,
         y1=plot.adfert$Y1,
         x0=plot.adfert$X,
         x1=plot.adfert$X,
         col=plot.adfert$col,
         lwd=1,
         lty=1)
points(x=plot.adfert$X,
       y=plot.adfert$Estimate,
       col=plot.adfert$col,
       pch=plot.adfert$pch,
       bg=plot.adfert$bg,
       cex=1.5)
axis(side=1,at=plot.adfert[1:8,]$Index,
     labels=plot.adfert[1:8,]$Model,
     cex.axis=1.25, tick=TRUE)




plot(1,
     pch="",
     ylim=c(1.15*min(plot.contraception$Y0),1.15*max(plot.contraception$Y1)), 
     xlim=c(.5,max(plot.contraception$Index)+.5), 
     xaxt="n",
     xlab="",
     ylab="Coefficient Estimate",
     main="Contraceptive Prevalence",
     cex.main=1.25,
     adj=1)
abline(h=0,lty="dotdash")
abline(v=3.5,lty="solid")

segments(y0=plot.contraception$Y0,
         y1=plot.contraception$Y1,
         x0=plot.contraception$X,
         x1=plot.contraception$X,
         col=plot.contraception$col,
         lwd=1,
         lty=1)
points(x=plot.contraception$X,
       y=plot.contraception$Estimate,
       col=plot.contraception$col,
       pch=plot.contraception$pch,
       bg=plot.contraception$bg,
       cex=1.5)
axis(side=1,at=plot.contraception[1:8,]$Index,
     labels=plot.contraception[1:8,]$Model,
     cex.axis=1.25, tick=TRUE)


plot(1,
     pch="",
     ylim=c(1.15*min(plot.GII$Y0),1.15*max(plot.GII$Y1)), 
     xlim=c(.5,max(plot.GII$Index)+.5), 
     xaxt="n",
     xlab="",
     ylab="Coefficient Estimate",
     main="UN Gender Inequality Index",
     cex.main=1.25,
     adj=1)
abline(h=0,lty="dotdash")
abline(v=3.5,lty="solid")

segments(y0=plot.GII$Y0,
         y1=plot.GII$Y1,
         x0=plot.GII$X,
         x1=plot.GII$X,
         col=plot.GII$col,
         lwd=1,
         lty=1)
points(x=plot.GII$X,
       y=plot.GII$Estimate,
       col=plot.GII$col,
       pch=plot.GII$pch,
       bg=plot.GII$bg,
       cex=1.5)
axis(side=1,at=plot.GII[1:8,]$Index,
     labels=plot.GII[1:8,]$Model,
     cex.axis=1.25, tick=TRUE)



plot(1,
     pch="",
     ylim=c(1.15*min(plot.alcohol$Y0),1.15*max(plot.alcohol$Y1)), 
     xlim=c(.5,max(plot.alcohol$Index)+.5), 
     xaxt="n",
     xlab="",
     ylab="Coefficient Estimate",
     main="Alcohol Consumption",
     cex.main=1.25,
     adj=1)
abline(h=0,lty="dotdash")
abline(v=3.5,lty="solid")

segments(y0=plot.alcohol$Y0,
         y1=plot.alcohol$Y1,
         x0=plot.alcohol$X,
         x1=plot.alcohol$X,
         col=plot.alcohol$col,
         lwd=1,
         lty=1)
points(x=plot.alcohol$X,
       y=plot.alcohol$Estimate,
       col=plot.alcohol$col,
       pch=plot.alcohol$pch,
       bg=plot.alcohol$bg,
       cex=1.5)
axis(side=1,at=plot.alcohol[1:8,]$Index,
     labels=plot.alcohol[1:8,]$Model,
     cex.axis=1.25, tick=TRUE)



plot(1,
     pch="",
     ylim=c(1.15*min(plot.malesmoking$Y0),1.15*max(plot.malesmoking$Y1)), 
     xlim=c(.5,max(plot.malesmoking$Index)+.5), 
     xaxt="n",
     xlab="",
     ylab="Coefficient Estimate",
     main="Tobacco Consumption (Male)",
     cex.main=1.25,
     adj=1)
abline(h=0,lty="dotdash")

segments(y0=plot.malesmoking$Y0,
         y1=plot.malesmoking$Y1,
         x0=plot.malesmoking$X,
         x1=plot.malesmoking$X,
         col=plot.malesmoking$col,
         lwd=1,
         lty=1)
points(x=plot.malesmoking$X,
       y=plot.malesmoking$Estimate,
       col=plot.malesmoking$col,
       pch=plot.malesmoking$pch,
       bg=plot.malesmoking$bg,
       cex=1.5)
axis(side=1,at=plot.malesmoking[1:8,]$Index,
     labels=plot.malesmoking[1:8,]$Model,
     cex.axis=1.25, tick=TRUE)




plot(1,
     pch="",
     ylim=c(1.15*min(plot.femsmoking$Y0),1.15*max(plot.femsmoking$Y1)), 
     xlim=c(.5,max(plot.femsmoking$Index)+.5), 
     xaxt="n",
     xlab="",
     ylab="Coefficient Estimate",
     main="Tobacco Consumption (Female)",
     cex.main=1.25,
     adj=1)
abline(h=0,lty="dotdash")

segments(y0=plot.femsmoking$Y0,
         y1=plot.femsmoking$Y1,
         x0=plot.femsmoking$X,
         x1=plot.femsmoking$X,
         col=plot.femsmoking$col,
         lwd=1,
         lty=1)
points(x=plot.femsmoking$X,
       y=plot.femsmoking$Estimate,
       col=plot.femsmoking$col,
       pch=plot.femsmoking$pch,
       bg=plot.femsmoking$bg,
       cex=1.5)
axis(side=1,at=plot.femsmoking[1:8,]$Index,
     labels=plot.femsmoking[1:8,]$Model,
     cex.axis=1.25, tick=TRUE)






# Title for whole thing
mtext("Religious Stratification Tests (Circles = Religious, Diamonds = Secular)", side=3, line=1, outer=TRUE, cex=1.15, font=2)

dev.off()

